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F33615-96-C-5070  SBIR  Topic  AF96-150:  Phase  I 


Final  Report 


Summary  Report 

WL/MLBM  solicited  Phase  I  SBIR  proposals  to  develop  a  capability  for  3D 
boundary  element  analysis  for  composite  joints  with  discrete  damage. 
Fracture  Analysis  Consultants,  Inc.  (FAC),  with  its  extensive  experience 
with  boimdary  element  technology  and  existing  boimdary  element 
softwEire,  responded  to  the  solicitation,  and  was  granted  a  six  month 
contract  to  demonstrate  the  feasibility  of  such  a  system. 

Our  overall  goal  for  was  to  demonstrate  expertise  in  BEM  methods 
appropriate  to  the  detailed  stress  analysis  of  arbitrary  bolted  or  bonded 
composite  joints  containing  cracks  and  contact  siu^aces.  There  were  three 
tasks  in  the  Phase  I  workstatement.  All  three  have  been  completed; 

Task  1.  We  augmented  the  Galerkin-based,  2D,  orthotropic  BEM  program 
developed  for  the  original  program  solicitation  to  include 
hypersingular  BEM  traction  equations  so  that  a  symmetric 
coefficient  matrix  could  be  formulated.  We  believe  that  this 
orthotropic  S3mimetric-Galerkin  code  is  unique. 

Task  2.  Using  the  code  developed  in  Task  1,  we  solved  an  example 

laminate  problem  [Pagano  1978],  to  show  that  the  symmetric- 
Galerkin  (SG)  approach  can  3deld  accurate  interlaminar  stress  . 
predictions  near  free  surfaces. 

Task  3.  We  demonstrated  a  capability  to  perform  3D  contact  analysis 

using  our  existing  3D  BEM  program.  The  approach  is  a  penalty 
function  technique  derived  from  a  current  capabilities  for  fluid 
driven  fracture,  and  can  eiIso  be  used  for  the  proposed  Phase  II  3D 
simulation  system. 

In  addition  to  these  tasks,  considerable  resources  were  spent  on  three 
additional  tasks  that  further  our  goal  of  demonstrating  expertise,  and  are 
preparatory  for  a  potential  Phase  II  effort. 

Task  4.  We  have  added  fracture  analysis  capabilities  to  the  2D  SG  code  of 
Task  1.  We  believe  that  this  is  the  first  SG  fracture  formulation 
for  elasto-statics,  and  certainly  the  first  for  orthotropic  materials. 

Task  5.  We  have  developed  software  to  evaluate  the  3D  anisotropic 

Green's  function.  This  is  a  generalization  of  the  work  of  Sales  and 
Gray  [1996].  In  this  approach  residue  calculus  is  used  to 
transform  the  Green's  function  to  a  ratio  of  polynomials.  A 
symbolic  algebra  program  is  used  to  generate  polynomial 
coefficients.  Newton's  method  is  used  to  find  the  roots  of  the 
polynomial.  Our  implementation  can  handle  multiple  orthotropic 
materials;  no  precomputations  are  required. 
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Task  6.  We  have  developed  a  fully  3D  Galerkin  BEM  progr£im.  Coding  for 
this  program  has  been  completed,  we  are  currently  in  the  testing 
and  debugging  phase.  This  code  includes  linear  and  quadratic 
order  elements,  and  has  been  interfaced  with  our  FRANC3D 
program  that  we  use  to  generate  models,  meshes,  emd  boimdary 
conditions. 

In  summary,  we  have  successfully  completed  the  work  tasks  specified  in 
our  propos^.  In  addition,  we  have  gone  substantially  beyond  the  original 
Phase  I  objectives  towards  creating  a  capability  for  state-of-the-art  3D 
boundary  element  analysis  of  composite  joints  with  discrete  damage. 


1.  Introduction  and  Background 

WL/MLBM  is  developing  improved  capability  for  detailed  stress  analysis  of 
composite  bolted  and  bonded  joints.  Such  joints  are  crucial  to  the  success 
of  composite  structiires  because  they  are  sites  for  initiation  of  a  variety  of 
damage  t3q)es,  and  contribute  substantially  to  weight  and  to  design  and 
manufactiiring  resources. 

Geometry  changes,  ply  drops,  and  free-edges  common  to  such  joints  create 
stress  concentrations  [Pagano,  1978]  and  theoretical  stress  singularities 
[Wang  &  Crossman,  1977]  which  render  simple  lamination  theory 
inapplicable.  A  mathematical  model  whose  characteristic  dimension  is 
lamina  thickness  is  required  to  approach  these  problems  [Chamis,  1989]. 
Although  there  have  been  some  successful  solution  techniques  based  on 
more  rigorous  analytical  approaches  at  this  scale,  these  are  restricted  in 
applicability  because  they  lack  the  geometrical  and  boundary  condition 
generality  of  numerical  methods  such  as  the  finite  element  method  (FEM). 

However,  there  has  been  a  historical  difficulty  in  using  the  FEM  for 
detailed  stress  analysis  of  laminates.  This  is  because  of  the  intrinsic 
difficulty  with  the  FEM  of  simultaneously  satisfying  displacement  and 
traction  compatibility  across  lamina  interfaces  while  also  providing  the 
strain  discontinuity  as  required  by  the  orthotropic  constitutive  laws  of  the 
lamina  materials  [Bogdanovich,  1992]. 

The  boundary  element  method  (BEM)  appears  to  be  an  attractive 
alternative  numerical  approach  for  this  application  because  it 
automatically  satisfies  displacement  and  traction  compatibility  across 
lamina  interfaces  while  naturally  reproducing  the  required  strain 
discontinuity. 

WL/MLBM  solicited  SBIR  proposals  for  the  development  of  a  3D  BEM 
capability  both  to  compare  with  in-house  variational  methods  of  solution, 
and  to  improve  efficiency  and  applicability.  Fracture  Analysis 
Consultants,  Inc.  (FAC),  with  its  extensive  experience  with  boxmdary 
element  teclmology  and  existing  boundary  element  software,  responded  to 
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the  solicitation,  and  was  granted  a  six  month  Phase  I  contract  to 
demonstrate  the  feasibility  of  such  a  system.  This  is  the  final  report  for 
this  Phase  I  effort. 


1.1  Phase  I  Technical  Objectives 

The  overall  goal  of  Phase  I  is  to  demonstrate  expertise  in  BEM  methods 
appropriate  to  the  detailed  stress  analysis  of  arbitrary  bolted  or  bonded 
composite  joints  containing  cracks  and  contact  simfaces. 

Technical  objectives  outlined  in  the  proposal  were; 

Objective  1:  Modification  of  a  new  Galerkin  2D  BEM  code  to  generalized 
plane  strain  and  to  symmetric  formulation  so  that  we  can  solve  the 
specified  test  problem  referenced  in  the  solicitation  [Pagano,  1978].  This 
will  demonstrate  the  feasibility  of  the  symmetric  Galerkin  boundary 
element  formulation. 

Objective  2:  Produce  a  fully  3D  solution  to  the  test  problem.  This  would 
then  enable  comparisons  for  accuracy  and  efficiency  to  the  predictions 
produced  from  Objective  1. 

Objective  3:  Modify  FRANC3D  (our  current  boundary  element  analysis 
system)  to  enable  solution  of  contact  stress  problems.  This  will  allow  us  to 
analyze  arbitrary  bolted  or  bonded  composite  joints  containing  contact 
surfaces. 

Objectives  1  and  3  have  been  fully  realized.  Objective  2  has  been  only 
partially  fulfilled.  Our  original  intention  for  Objective  2  was  to  modify  our 
existing,  collocation-based  finite  element  code  to  include  anisotropic 
material  capabilities.  However,  during  the  course  of  the  project  we  decided 
that  it  would  be  more  useful  in  the  long  run  to  develop  a  new  3D  Galerkin 
formulation  boimdary  element  code  (see  the  next  three  sections  for 
definitions  of  these  terms).  We  have  developed  this  new  capability,  but 
because  this  effort  was  more  than  that  of  the  original  obective  2,  we  have 
not  yet  been  able  to  perform  a  fully  3D  solution  of  the  test  problem. 

A  detailed  description  of  the  results  of  all  tasks  associated  with  these 
objectives  is  provided  in  Section  2.  As  noted  there,  even  though  objective  2 
had  been  only  partialy  fullfilled,  all  tasks  have  been  completed,  along  with 
a  number  of  additional  tasks  not  in  the  original  workstatement.  First, 
however,  we  provide  a  brief  backgroxmd  on  boimdary  element  formulations 
and  computational  fracture  analysis  therewith.  These  sections  proved  a 
technical  context  for  the  remainder  of  the  report. 


4 


1^  Boundary  Element  Formulations 


Conventionally,  the  boundary  element  is  formulated  as  a  collocation 
procedure.  That  is,  the  governing  equation 

“»(/’)  +  J  «i,(e)  dQ  =  J  U^^T,(Q)dQ  (1-1) 

is  forced  to  be  satisfied  at  some  number  of  nodal  points,  on  the  botmdary. 
The  equation  is  not  necessarily  satisfied  away  from  these  points.  While  P 
is  a  point,  the  integration  is  performed  over  a  surface  element  Q.  This  is 
shown  pictorially  in  Figure  1-la.  The  integration  is  relatively 
straightforward,  except  for  the  cases  where  the  P  node  is  adjacent  to  the  Q 
element.  As  the  traction  and  displacement  kernels,  U  and  T,  contain 
factors  of  l/r  and  respectively  (q*  and p*  being  the  Cartesian 

coordinates  of  the  source  and  field  points),  these  integrations  become 
singular,  and  special  care  must  be  taken  in  their  evaluation.  Techniques 
for  this  are  well  known,  and  can  be  found  in  any  number  of  standard  BEM 
textbooks  [e.g.,  Brebbia  and  Dominguez,  1989]. 

While  the  collocation  technique  forces  the  governing  equations  to  be 
satisfied  at  a  finite  number  of  points,  the  finite  element  technique  has  been 
very  successful  using  the  Galerkin  form  of  the  method  of  weighted 
residuals.  In  this  approach,  the  governing  equations  are  not  satisfied  at 
individual  points,  but  rather  in  an  average  sense  over  an  element's 
domain.  This  approach  can  be  applied  to  boundary  elements  [Maier  and 
Polizzotto  1987]  but  requires  two  integrations,  one  each  over  the  P  and  Q 
elements. 

Symbolically,  the  Galerkin  governing  equations  are 

J  xir,{P)u^{P)dP  +  f  iif,iP)\  T^jjiP,Q)  u^iQ)  dQdP 

,,  (1-2) 
=  ]V,iP)\u^jjT,{Q)dQdP 

where  y  ^^e  the  shape  functions  of  the  P  element.  This  is  shown  in  Figure 
1-lb. 


There  are  at  least  five  advantages  of  the  Galerkin  formulation: 

1.  The  technique  to  yields  much  more  accurate  internal  stresses  and 
displacements  close  to  the  boundaries  than  does  the  collocation 
approach.  Such  stresses  and  displacements  are  notoriously  bad  with 
collocation  giving  rise  to  the  so-called  "boundary  effect".  In  essence, 
since  the  governing  equations  are  satisfied  only  at  nodal  points, 
evaluating  stresses  or  displacements  at  points  near  the  boimdary  can 
give  rise  to  significant  errors.  The  recommendation  is  to  not  trust 
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values  evaluated  closer  to  the  boundary  than  the  local  characteristic 
nodal  spacing. 


Figure  1-1.  a)  normal  collocation  boundary  element  integration,  b) 
Galerkin  boundary  element  integration. 

Such  restrictions  on  evaluating  stresses  and  displacements  are  quite 
severe  in  the  present  case  where  we  are  concerned  with  inter-laminar 
stresses,  and  stress  fields  in  near  crack-tip  regions.  Therefore,  the 
Galerkin  formulation  seems  superior  for  these  cases. 

2.  The  technique  gives  one  the  ability  to  rigorously  treat  hypersingular 
equations  using  standard  continuous  elements,  collocation  requiring  a 
differentiable  approximation  (this  is  discussed  further  below). 

3.  A  more  accurate  treatment  of  comers,  especially  for  corners  having 
displacement  boimdary  conditions  and  corners  having  no  boundary 
conditions.  This  is  important  for  bi-material  interface  problems  and 
potentially  for  coupling  BEM  domains  to  FEM  domains. 

4.  In  conjimction  with  the  hypersingular  boundary  element  equations,  a 
symmetric  coefficient  matrix  can  be  formulated,  giving  rise  to 
significant  decreases  in  necessary  storage  and  processing  time. 

5.  A  sjnnmetric  formulation  give  the  possibility  to  couple  boundary 
elements  with  finite  elements  and  still  maintain  the  symmetric  FEM 
coefficient,  matrix. 


1^  Boundary  Elements  and  Fracture  Mechanics. 

A  naive  application  of  boundary  element  techniques  to  fracture  mechanics 
fails.  This  is  illustrated  in  Figure  l-2a.  Nodes  i  and  j,  while  being  on 
distinct  surfaces,  are  geometrically  coincident.  This  means  that  the 
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boundary  element  equations  for  these  two  points  will  be  identical,  leading 
to  a  rank  deficient  coefficient  matrix. 


The  remedy  for  this  situation  has  been  to  break  the  object  into  two 
separate  domains  by  way  of  a  fictitious  boimdary  [Rizzo  and  Sbippy  1968, 
Cruse  1969,  Blanford,  Ingrafifea,  and  Liggett  1981].  Separate  boundary 
element  coefficient  matrices  are  formulated  for  the  two  domains,  and 
equilibrium  and  compatibility  is  enforced  along  the  boundary.  This  is 
illustrated  in  Figure  l-2b. 


domain  I 

fictitious  boundaries 

node  j  ^ 

crack 

,  .  __  crack 

„  ,  domain  II 

(a) 

(b) 

Figure  1-2.  a)  a  naive  application  of  boundary  elements  to  a  fracture 
mechanics  problem,  b)  the  multi-domain  approach. 

The  multi-domain  approach  has  been  shown  to  be  unsatisfactory  for  at 
least  two  reasons.  First,  while  for  two  dimensional  problems  defining  a 
suitable  fictitious  boimdary  is  usually  straightforward,  doing  the  same  for 
an  arbitrary  three-dimensional  configuration  can  be  quite  difficult. 

Second,  this  approach  adds  additional  surfaces  in  the  crack  tip  region.  As 
mentioned  above,  while  the  governing  equations  are  satisfied  exactly  in 
the  interior  of  the  domain,  they  are  only  approximately  satisfied  near  the 
boundaries.  The  multi-domain  approach  therefore  places  additional 
approximations  in  the  crack-tip  region,  where  stress  gradients  are 
expected  to  be  the  highest,  and  where  usually,  the  most  accurate  results 
are  desired. 

Gray  and  his  associates  have  taken  a  different  approach  to  the  crack 
problem,  where  the  redundant  crack  face  equations  are  replaced  with 
equations  for  the  known  (usually  zero)  stress  conditions  on  the  crack  faces 
[Gray  1989,  Gray,  Martha,  and  Ingraffea  1990,  Lutz,  Gray,  and  Ingrafifea 
1991,  Martha  et.  al.  1992,  Gray  1993]. 

In  his  approach,  Eq.  (1-1)  is  differentiated  with  respect  to  P  to  yield  the 
following  equation  for  the  tractions, 

W  +  J  S^i(.P,Q)Uji(Q)dQ  =  I  W^,{P,Q)T,(,Q)dQ  (1-3) 

which  can  be  substituted  for  the  redundant  crack-face  equations. 


The  S  and  W  kernels  of  Eq.  (1-3)  contain  first  and  second  derivatives  of  the 
Green's  function,  and  thus  lead  to  so-called  hypersingular  equations. 

Gray  has  shown  how  these  can  be  evaluated  accurately  by  decomposing 
the  integrals  into  the  non-singular  parts,  integrated  numerically,  and 
simple  singular  parts,  evaluated  analjrtically. 

The  field  of  hypersingular  boundary  integral  equations  has  opened 
possibilities  for  new  solutions  to  various  problems  by  the  BEM  [see,  for 
example,  the  review  article  by  Krishnasamy,  Rizzo,  and  Rudolphi  1992]. 
There  are  a  number  of  hypersingular  approaches,  and  each  approach  has 
been  developed  independently  by  sever^  groups:  e.g.  Bonnet  and  Bui 
[1993],  Crouch  and  Selcuk  [1992],  Cruse  [1988],  Chang  and  Mear  [1995], 
Gray,  Martha,  and  Ingraffea  [1990],  Guimaraes  and  Telles  [1994],  Hong 
and  Chen  [1988],  and  Paulino  [1995].  Although  successful,  the  common 
problem  for  these  methods  is  that,  when  employed  in  conjunction  with  a 
collocation  approximation,  a  continuous  smoothness  constraint  is 
necessarily  imposed  on  the  boundary  in  the  neighborhood  of  the 
integration  point.  This  condition  adds  considerable  complexity  and/or  cost 
to  fractxire  analyses.  Historically,  we  have  satisfied  this  condition  by 
moving  the  collocation  points  from  the  element  nodal  points  to  points  on 
the  interior  of  the  elements,  with  the  corresponding  cost  of  significantly 
larger  system  of  equations. 

However,  if  Eq.  (1-3)  is  evaluated  with  a  Galerkin  formulation,  the 
resulting  equation, 

J  (P)  (/>)  <iP  +  f  y.  (/>) 1 (/>,  0  (0  dSdP 

,  ,  (1-4) 

=  J^(P)Jw'„(P,0r/0daiP 

can  be  evaluated  with  only  C°  continuity  between  elements. 


1.4  A  Symmetric  Galerkin  Formulation 

The  key  to  obtaining  a  sjmimetric  coefficient  matrix  are  the  properties  of 
the  kernel  functions  [Sirtori  et.  al.  1992], 

U(P,Q)  =  UiQ,P) 
nP,Q)  =  -T{Q,P) 

W{P,Q)  =  W{Q,P)  ^  ' 

SiP,Q)  =  -S(Q,P) 

combined  with  a  Galerkin  approximation,  and  appropriate  selection  of 
equations  for  each  degree  of  freedom.  The  displacement  equation  (1-2)  is 
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used  for  all  degrees  of  freedom  where  displacements  are  specified,  and  the 
traction  equation  (1-4)  is  used  for  the  remaining  degrees  of  freedom  where 
tractions  are  specified.  After  discretization,  the  resulting  set  of  equations 
can  be  written  in  the  following  block  matrix  notation. 


(1-6) 


where  the  k  and  u  subscripts  represent,  respectively,  the  known  and 
unknown  components  of  the  displacement  and  traction  vectors. 

Rearranging  Eq.  1-6  into  the  form  [A]{x}  =  {b},  and  multiplying  the  second 
row  (hypersingular  equations)  by  -1,  one  obtains. 


-//„«*  +  GijTJ 

,  ^21%  —  ^22^t  J 


(1-7) 


The  symmetry  of  the  coefficient  matrix,  G,,  =  G,^,  H22  =  HI2,  and  H^2  =  Gl^, 
now  follows  from  the  properties  of  the  kernel  fimctions,  Eq.  (1-5). 


2.  Workstatement  Tasks 

There  were  three  tasks  in  the  proposed  work  statement  of  this  project.  The 
first  two  were  aimed  at  demonstrating  that  a  symmetric  Galerkin 
boundary  element  formulation  would  work,  and  are  continuations  of  the 
analysis  included  in  our  Phase  I  proposal.  The  third  task  was  looking 
forward  to  potential  Phase  II  work,  and  was  to  investigate  the  possibility  of 
incorporating  contact  analysis  into  a  3D  boundary  element  analysis.  All 
three  tasks  have  been  completed. 


2.1  Task  I  -  Creation  of  a  2D  Symmetric  Galerkin  Code 

The  program  solicitation  for  this  project  requested  that  all  proposals 
contain  a  plane  strain,  thermal-mechanical  analysis  of  a  [0790°]s 
laminate  subjected  to  a  uniform  temperature  change  of  -100°F.  At  the 
time  of  the  solicitation,  we  had  no  two-dimensional  BEM  capabilities 
(while  we  did  have  three-dimensional  capabilities).  Therefore,  we 
developed  a  new  two-dimensional  code  for  the  purpose  of  responding  to  the 
solicitation. 

This  new  code  was  a  Galerkin  formulation,  but  not  symmetric.  That  is,  it 
used  the  displacement  equation  (1-2)  but  did  not  incorporate  the  traction 
equation  (1-4),  We  proposed,  as  Task  I,  to  augment  our  2D  code  to 
incorporate  a  sjmimetric  formulation,  and  to  include  capabilities  for 
generalized  plane  strain. 
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A  symmetric  formulation  has  been  implemented.  There  is  essentially  no 
difference  in  computed  results  between  the  symmetric  and  unsymmetric 
results.  The  time  to  formulate  elements  is  essentially  the  same  for  the  two 
programs,  while  the  solution  time  for  the  S3nnmetric  case  is  approximately 
one  half  that  of  the  unsymmetric  version. 

The  generalized  plane  strain  capability  is  a  near  trivial  extension  of  the 
way  we  implemented  thermo-mechanical  loading.  Assume  that  the  plane 
of  analysis  is  the  y-z  plane.  Given  the  full  three  dimensional  constitutive 
relationships, 
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and  the  initial  strain  vector  due  to  a  constant  temperature  change, 

[a,AT  a^AT  a^AT  0  0  (2-2) 


the  expression  for  the  in-plane  stresses  is 
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a^AT .  (2-3) 


In  our  code,  the  body  stress  terms 


(2-4) 


are  added  to  accoimt  for  the  thermal  strains.  A  similar  approach  can  be 
taken  for  generalized  plane  strain.  However,  in  this  case  the  body  stress 
terms  are 


which  is  numerically  the  same  expression  as  2-4  with  =  a,  =  0 ,  and 

e,  =  . 


10 


F33615-96-C-5070  SBIR  Topic  AF96-150:  Phase  I 


Final  Report 


2^  Task  n  -  Solution  of  a  Test  Problem 

To  demonstrate  the  capabilities  added  in  Task  I,  a  problem  originally 
examined  by  Pagano  [1978]  was  reanalyzed  as  Task  II.  The  problem  is  the 
[0°/90°ls  laminate  subjected  to  a  uniform  x  strain,  as  shown  in  Figure  2-1. 


Figure  2-1.  A  [0790°]s  laminate  subjected  to  uniform  strain  in  the  x  - 

direction. 

The  laminate  is  modeled  using  generalized  plane  strain,  and  only  one 
quarter  of  the  cross-section  was  modeled,  as  shown  in  Figure  2-2. 


£;i=20xl0*psi 

£^2  =  F33  =  2.1xl0®psi 

Gi2  =  G"i3  =  G23  =  0.85xl0^psi 

Vi2  =  V,3  =  V23  =  0.21 


Figure  2-2.  The  quarter  symmetry  model  and  material  properties  for  the 

test  problem. 


The  computed  inter-lamina  stress  distributions  along  the  0790°  interface 
are  shown  in  Figures  2-3  and  2-4  (these  are  directly  comparable  to  Figures 
7  emd  8  of  the  Pagano  paper). 

As  cem  be  seen  in  Figure  2-3,  there  is  very  good  correlation  between  our 
computed  normal  stress  distribution  and  that  reported  by  Pagano  (note, 
the  Pagano  results  plotted  here  were  determined  graphically  from  the 
figures  in  the  paper  and  thus  are  only  approximate). 
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The  comparison  for  the  shear  stresses  is  not  quite  as  good.  As  shown  in 
Figure  2-4,  the  present  analysis  predicts  higher  peak  shear  stresses  than 
the  reference  analysis.  We  have  not  investigated  the  effect  of  mesh 
refinement. 


Figure  2-3.  Comparison  of  the  normal  inter-lamina  stresses  for  the 
present  analysis  and  that  of  Pagano  [1978]. 


y/b 

Figure  2-4.  Comparison  of  the  shear  inter-lamina  stresses  for  the  present 

analysis  and  that  of  Pagano  [1978]. 
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2.3  Task  m  -  3D  BEM  Contact  Analyses 

An  ultimate  objective  of  this  project  is  to  develop  a  detailed  analysis 
capability  for  mechanically  fastened  joints  in  composite  materials.  In 
such  a  joint,  bearing  stresses  from  bolts  play  a  major  role  in  the  load 
transfer  through  the  joint,  and  must  be  modeled  properly  in  a  detailed 
local  stress  analysis.  Such  an  analysis  must  include  elastic  deformation  of 
the  bolt  and  evaluation  of  contact  zones.  In  the  Phase  I  effort,  we  have 
demonstrated  this  capability  as  a  straightforward  technology  transfer 
from  our  current  coupled  fluid/structure  interaction  [Sousa  et  al,  1993] 
and  non-linear  cohesive  fracture  capabilities  [Bittencourt  and  Ingraffea, 
1995]. 

An  important  application  of  our  current  boundary  element  program  is  the 
analysis  of  fluid  driven  fracture  (hydraulic  fracture).  In  this  case, 
cracking  is  driven  by  pressurized,  non-Newtonian  fluids  inside  the  crack. 
This  is  a  highly  nonlinear,  coupled  fluid/structure  interaction  problem. 

The  nonlinearity  arises  because  the  flow  in  the  crack  is  a  cubic  function  of 
the  crack  opening.  The  solid,  however,  remains  linear  and  elastic  in 
behavior. 

We  solve  this  problem  by  using  the  boundary  element  program  to  generate 
a  flexibility  matrix  for  the  crack  surface,  which  relates  the  displacement 
at  any  point  of  the  crack  surface  to  a  load  at  any  other  point  on  the  surface. 
In  essence,  this  is  using  the  boundary  element  program  to  create  a  super 
element.  Once  the  flexibility  matrix  has  been  generated,  we  use  a 
nonlinear  solution  algorithm  for  solving  the  fluid-structure  interaction 
equations.  The  important  point  here  is  that,  at  each  iteration  in  the 
solution,  we  only  need  to  perform  a  matrix  multiply  to  incorporate  the 
effects  of  the  solid  in  the  solution. 

We  have  adapted  our  fluid/structure  interaction  capabilities, 
incorporating  a  penalty  function  approach,  for  solving  contact  problems. 
The  general  contact  problem  has  been  modeled  in  3D  using  FRANC3D 
[Wawrzynek  1991]  and  our  existing  BEM  program,  BES.  The  contact 
surface  is  modeled  as  two  distinct  siufaces  that  share  the  same  geometric 
description.  These  surfaces  are  meshed  with  the  same  surface  mesh. 
Multiple  solutions  are  obtained  during  the  solution  process;  one 
corresponds  to  the  far  field  applied  loading,  while  the  others  correspond  to 
rniit  tractions  applied  to  each  node  pair  on  the  contact  surface.  This 
provides  a  set  of  influence  coefficients  that  can  be  used  to  determine  the 
equilibrium  displacements  and  tractions  due  to  unit  tractions  applied  to 
the  contact  surfaces. 

An  iterative  scheme  is  used  to  determine  the  equilibrium  tractions  on  the 
contact  surface.  The  algorithm  is  summarized  below; 
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1.  Assume  that  all  tractions  are  zero  on  the  contact  surface,  and  the 
amoimt  of  overlap  due  to  the  far  field  loading  is  determined; 

2.  At  the  node  where  maximum  overlap  occurs,  a  traction  is  applied.  The 
traction  is  computed  by  comparing  the  displacements  between  the  node 
pair  on  the  contact  surfaces  for  the  cases  of  applied  far  field  loading  and 
a  unit  traction  on  the  node  pair  (unit  tractions  are  applied  in  the 
opposite  sense  to  the  two  nodes); 

3.  Based  on  the  ratio  of  the  displacements  and  a  numerical  damping 
(releixation)  factor,  a  traction  value  is  computed  and  applied  to  the  node 
pair; 

4.  The  solution  for  the  far  field  loading  is  then  superimposed  with  the 
solution  for  the  applied  traction  on  the  contact  surface  assuming  linear 
elasticity; 

5.  The  total  displacement  and  traction  field  is  then  recomputed; 

6.  A  search  is  made  again  for  the  node  pair  with  the  maximum  overlap, 
and  the  process  is  repeated,  summing  up  all  solutions  for  all  non-zero 
tractions  and  computing  new  displacements  and  tractions  at  each 
iteration. 

As  an  example,  a  model  similar  to  that  analyzed  by  Yogeswaren  and  Reddy 
[1988]  was  considered.  It  should,  however,  be  emphasized  that  while 
Yogeswaren  and  Reddy  performed  plane  strain  calculations,  the  present 
aniysis  is  a  fully  three-dimensional  treatment.  They  were  modeling 
experiments  performed  by  Joh  [1986].  The  model  consists  of  an  isotropic 
alumimun  plate  with  a  hole  and  pin.  The  pin  is  fixed  and  tension  is 
applied  to  the  plate.  The  model,  material  parameters  and  loading  are 
shown  in  Figure  2-5.  Appropriate  restraints  are  applied  to  prevent  rigid 
body  motion.  The  plate  and  pin  are  0.06"  thick.  Friction  on  the  contact 
surface  is  not  modeled.  Yogeswaren  and  Reddy  modeled  a  clearance  fit  pin 
with  0.005"  between  the  pin  and  the  hole,  and  they  included  friction  on  the 
contact  surfaces;  both  a  static  and  d3niamic  coefficient  of  fHction  were 
used.  In  the  present  analysis,  a  slip-fit  was  assumed,  and  fictional  effects 
were  ignored. 

For  the  first  case  analyzed,  the  pin  has  the  same  material  properties  as  the 
plate.  The  displacements  due  to  far  field  loading  alone  are  shown  in  Figure 
2-6a  and  the  actual  displacements  after  accounting  for  contact  are  shown 
in  Figure  2-6b.  Normal  tractions  on  the  contact  surface  are  shown  in 
Figure  2-7.  Zero  degrees  corresponds  to  position  A  in  Figure  2-5.  The 
minimum  traction  (maximum  compression)  exists  at  about  65  degrees. 

This  is  edso  indicated  by  the  color  contour  plot  of  minimum  principal  stress 
(Figure  2-8). 
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Figiire  2-6.  (a)  displaced  shape  due  to  far  field  loading  only,  and  (b) 
displaced  shaped  after  accounting  for  contact. 
Magnification  factor  is  100. 
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Figure  2-8.  Color  contours  of  the  minimum  principal  stress  also  shows 
that  the  maximum  compressive  contact  stresses  are  at  about  65  degrees. 

Magnification  factor  is  100. 

The  minimum  value  of  normal  traction  are  substanitally  less  than  the 
experimental  experimental  values  reported  Joh  (about  -35  ksi).  It  is 
imcleeir  from  the  Yogeswaren  and  Reddy  paper  whether  the  reported  load 
value  of  1210  lbs.  was  the  total  load  applied  to  the  experiments,  or  the  total 
load  applied  to  their  finite  element  mesh,  which  models  one  half  of  the 
experiment.  In  retrospect,  it  appears  that  the  later  is  the  case,  while  we 
assumed  that  the  reported  load  was  the  total  experimental  load,  and  only 
applied  605  lbs.  to  our  half  symmetry  analysis.  Furthermore,  Rao  [1978] 
has  shown  that  the  amount  of  clearance  and  friction  can  significantly 
effect  the  contact  stress  distribution  and  location  of  separation.  Therefore, 
direct  comparisons  between  our  analysis  and  the  experiments  by  Joh 
cannot  be  made. 

A  second  model  was  analyzed  with  a  stifFer  "steel"  pin  (E=30,000  ksi).  The 
normal  traction  on  the  contact  surface  for  this  model  is  shown  in  Figure  2- 
9.  The  trend  of  the  curve  is  similar  to  the  aluminum  pin,  although  the 
minimum  is  shifted  to  about  60  degrees. 

In  both  cases,  the  minimum  normal  traction  on  the  contact  surface  is  not 
at  position  A  (0  degrees),  because  the  pin  is  close  to  the  fixed  end  of  the 
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plate.  A  third  model  was  created  by  extending  the  plate  3"  to  the  left, 
leaving  the  pin  much  further  from  the  fixed  end.  The  steel  pin  was  used  in 
this  model.  The  normal  traction  on  the  contact  surface  for  this  model 
(Figure  2-9)  shows  that  the  normal  traction  is  almost  constant  from  0  to 
about  60  degrees  after  which  it  decays  to  zero. 


Figure  2-9.  Contact  normal  stress  for  an  aluminum  pin,  steel  pin, 
and  a  steel  pin  with  an  extended  plate. 


3.  Additional  Tasks  not  in  the  Workstatement 

There  were  three  tasks  not  in  the  workstatement  of  our  Phase  I  proposal 
for  which  we  spent  a  considerable  amount  of  effort  during  the  course  of 
this  project.  In  fact,  these  efforts  represent  the  majority  of  our  effort. 

First,  we  extended  o\ir  two-dimension  symmetric-Galerkin  BEM  code  to 
include  single  domain  fracture  analysis,  for  orthotropic  materials.  Second, 
in  anticipation  of  three-dimensional  developments,  we  formulated  an 
efficient  technique  for  evaluating  3D  anisotropic  Green's  functions. 
Finally,  we  have  formulated  and  implemented  a  three-dimensional 
anisotropic  Galerkin  (but  not  yet  symmetric)  BEM  code. 


3.1  2D  Symmetric  Galerkin  Fracture  Mechanics 

A  symmetric  Galerkin  algorithm  for  fracture,  in  the  context  of  two- 
dimensional  potential  theory,  was  recently  presented  by  Gray, 
Balakrishna,  and  Kane  [1995].  This  method,  effectively  a  combination  of 
displacement  discontinuity  [Crouch  and  Starfield  1983]  and  the 
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hypersingular  BEM  [e.g.,  Gray  1990]  ideas,  employs  the  hypersingular 
equation  on  the  crack  surface  as  described  above  in  Section  1.3.  For  the 
present  project,  we  have  explored  techniques  to  improve  the  computational 
efficiency,  and  have  extended  the  approach  to  plane  orthotropic  elasto- 
statics.  Some  of  this  work  has  been  presented  in  Gray  and  Paulino  [1996]. 


3.1.1  Formulation 

Consider  a  2D,  linear-elastic,  orthotropic  body  which  contains  a  crack,  as 
illustrated  in  Rgure  3-1.  The  body  has  a  botmdary  F  =  r„  +  -f  +  F^, 

which  represents,  respectively,  the  portion  of  the  boundary  over  which 
displacements  are  specified,  tractions  are  specified,  and  the  two  sides  of 
the  crack  (over  which  the  tractions  are  assumed  to  be  known). 


n  ^ 

7 

/ 

r. 

/  I. 

f; 

r„ 

'777777. 

'7777 

Figure  3-1  A  generic  body  containing  a  crack. 


In  the  hypersingular  approach  to  solving  fracture  problems  described  in 
Section  1.3,  independent  equations  for  solving  for  the  displacements  on 
both  sides  of  the  crack  are  obtained  by  writing  both  standard  and 
h3q)ersingular  equations.  This  violates  the  symmetric-Galerkin 
procedure,  as  the  boundary  condition  on  the  crack  surfaces  are  assumed  to 
be  specified  tractions  and,  consequently,  only  the  traction  equation  should 
be  used.  Symmetry  of  the  coefficient  matrix  is  therefore  not  possible  with 
the  hypersingular  approach.  Nevertheless,  taking  this  as  a  starting  point, 
the  resxilting  system  of  equations  can  be  written  in  block-matrix  form. 
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The  blocking  strategy,  according  to  Eq.  (3-1),  is  as  follows:  the  first  row 
and  column  is  associated  with  the  outer,  non-crack  boundary  F„  h-Fj,  while 
subscripts  2  and  3  refer  to  the  two  sides  of  the  crack  F^  and  FJ", 
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respectively.  The  vector  of  unknowns  on  the  outer  boundary  is,  in  general, 
a  mixture  of  displacements  and  tractions,  and  is  denoted  by  .  The 
corresponding  vector  of  prescribed  boundary  values  is  indicated  by  'P,. 

The  second  and  third  rows  correspond  to  the  equations  for  traction  and 
displacement,  respectively,  written  on  the  crack  surface  (again,  it  assumed 
that  the  traction  is  specified  on  the  crack  surface). 


Note  that,  although  the  coefficient  matrix  in  Eq.  (3-1)  is  not  symmetric,  the 
upper  2x2  principal  submatrix  is,  i.e.  =  ///,,  =  Z/^,  and  =  ZZ/j. 

This  is  a  consequence  of  the  basic  symmetric-Galerkin  procedure,  the 
traction  being  the  appropriate  equation  choice  on  the  crack.  In  addition  to 
the  symmetry,  it  is  also  important  to  observe  that  Z/j,  =  -Z/jj  and  Z/jj  =  -ZZjj, 
a  consequence  of  the  change  in  simface  orientation  in  the  integration  over 
the  two  sides  of  the  crack.  It  can  be  shown  that  a  sjmametric  coefficient 
matrix  results  from  changing  variables  on  the  crack  from  displacement  to 
the  jump  in  displacement  (i.e.,  displacement  discontinuity). 


Am  =  Oj  -  Mj 

Zm  =  M2  +M3 


(3-2) 


With  this  transformation,  the  left  hand  side  of  Eq.  (3-1)  takes  the  form 


'Hn 

^,2 

0  ■ 

Qi 

^12 

0 

Am 

_Z/3. 

III 

Zm 

(3-3) 


where  I  is  the  identity  matrix.  It  therefore  suffices  to  solve  the  smaller 
symmetric  2x2  block  system  for  the  unknowns  If  necessary,  Zm, 

and  hence  the  actual  crack  face  displacements  can  be  calculated  in  a  post¬ 
processing  step. 


3.1.2  Example  Problems 

Two  examples  are  presented  to  validate  the  2D,  ortho  tropic,  symmetric- 
Galerkin  fracture  capability.  Each  of  these  examples  involves  a  single 
interior  crack,  one  aligned  with  the  material  axes  (Figure  3-2a)  and  one 
inclined  at  a  45°  to  the  axes  (Figure  3-2b).  For  both  examples,  plane  stress 
conditions  are  assumed.  The  crack  has  length  2a  =  0.4,  the  plate  has 
height  2H  =  2.0  and  width  2W  =  1.0.  A  remote  traction  of  10.0  MPa  was 
applied.  The  elastic  constants  used,  which  correspond  to  average 
homogenized  properties  for  fiberglass,  are  [Ghandi  1972]: 

E^  =48.26  GPa 
£2  =  17.24  GPa 
/i,2  =  6.89  GPa 
v,2=  0.291 
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The  boundary  element  results  for  these  examples  are  compared  to  finite 
element  solutions  from  the  program  FRANC2D  [Wawrzynek  1991,  Boone, 
Wawrzynek,  and  Ingraffeal987].  The  parameter  chosen  for  comparison 
purposes  is  the  displacement  discontinuity  Am  along  the  crack. 

10  MPa  10  MPa 


Figure  3-2.  The  two  example  problems  used  to  validate  the  2D,  symmetric- 
Galerkin  boundary  element  capability. 


The  finite  element  solutions  used  approximately  270  quadratic  order 
elements  with  about  1560  degrees-of-freedom.  Quarter-point  elements  are 
used  to  model  the  stress  singularity  at  the  crack  tips. 

For  the  boimdary  element  solutions,  the  outer  boundary  of  the  plate  has 
been  uniformly  discretized  with  60  standard  linear  elements.  Each  face  of 
the  crack  has  been  discretized  with  8  linear  boundary  elements,  with  some 
grading  towards  finer  elements  towards  the  crack  tip.  Such 
discretization,  especially  of  the  crack  faces,  is  much  too  coarse  to  expect 
highly  accimate  results,  especially  in  the  crack  tip  region  where 
displacement  gradients  are  high.  However,  it  was  felt  that  if  meaningful 
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comparisons  to  the  FEM  results  could  be  made,  this  would  be  sufficient  to 
validate  the  approach. 

Table  3-1  shows  the  results  for  the  crack  opening  displacement  (COD)  for 
the  crack  of  example  1.  Because  of  symmetry,  only  the  results  for  the  left 
portion  of  the  crack  are  computed.  The  COD's  obtained  by  both  the  BEM 
and  FEM  analyses  Eire  in  reasonably  good  agreement,  the  largest 
discrepancy  being  less  than  5%.  For  this  problem,  the  crack  sliding 
displacement  (CSD)  should  be  zero.  The  numerical  values  for  both  BEM 
and  FEM  are  of  the  order  10  ®  or  less,  indicating  that  consistent  solutions 
have  been  obtained. 

Table  3-1.  Crack  opening  displacement  for  Example  1.  The  solution  is  for 
one  half  of  the  crack  starting  at  the  left  crack  tip. 


y 

SG-BEM 
Am^  (xlO-^) 

FEM 

Am,  (xlO-®) 

difference  | 
%  1 

1.000 

0.000000 

0.000000 

1.000 

0.227131 

0.230606 

0.500 

1.000 

0.391940 

0.411062 

0.450 

1.000 

0.439092 

0.460448 

0.500 

1.000 

0.453433 

0.476836 

4.91 

Table  3-2  shows  the  results  for  the  displacement  discontinuity  for  the 
entire  lenght  of  the  crack  of  example  2.  As  before,  the  results  for  Am 
obtained  by  both  the  BEM  and  FEM  are  in  reasonable  agreement,  the 
largest  discrepancy  being  about  5%.  However,  the  results  for  Au^  show 
much  larger  relative  discrepancies,  especially  for  the  nodes  near  the  crack 
tips.  It  is  quite  likely  that  for  this  case  neither  the  BEM  or  FEM  solutions 
are  converged.  However,  for  the  present  purposes,  these  analyses 
demonstrate  that  our  extensions  of  the  symmetric  Galerkin  boundary 
element  formulation  for  fracture  problems  is  valid. 

Table  3-2.  The  displacement  discontinuities  for  Example  2. 


X 

■ 

BEM 

FEM 

Am,(%) 

>  1 
s  U 

Am,(x10-') 

Am,(x10-') 

Am,(x10"') 

Am,(x10-') 

0.359 

0.859 

0.000000 

0.000000 

0.000000 

0.000000 

0.00 

0.00 

0.376 

0.876 

-.049903 

0.163476 

-.065051 

0.164245 

23.29 

0.47 

0.429 

0.929 

-.105063 

0.283705 

-.113276 

0.298861 

7.25 

5.07 

0.465 

0.965 

-.125796 

0.318479 

-.139346 

0.334257 

9.72 

4.72 

0.500 

1.000 

-.132743 

0.329145 

-.145958 

0.345614 

9.05 

4.77 

0.535 

1.035 

-.125730 

0.318613 

-.139309 

0.334217 

9.75 

4.67 

0.571 

1.071 

-.104960 

0.283939 

-.113400 

0.298794 

7.44 

4.97  1 
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0.624 

1.124 

-.049838 

0.163708 

-.064930 

0.164171 

23.24 

0.28 

0.641 

1.141 

0.000000 

0.000000 

0.000000 

0.000000 

0.00 

0.00 

3^  Effective  Evaluation  of  the  3D  Anisotropic  Green's  Function 

3.2.1  Formulation 

A  necessary  component  of  a  practical  program  for  performing  3D  BEM 
stress  analysis  for  orthotropic  materials  is  a  routine  for  efficiently 
evaluating  the  anisotropic  Green's  function.  This  function,  which  can  be 
stated  as: 


where  the  Christoffel  matrix,  K,  is  defined  in  terms  of  the  elastic 
constants,  C,  of  the  material, 

and,  r  =  l|r-yl|.  (3-5) 

l,m 

The  integration  path  is  the  unit  circle  in  the  plane  that  is  perpendicular  to 
the  vector  x-y, 


=  S\x,y)  =  {^ € 51'  IH  =  l,^.(x-y)  =  O}.  (3-6) 

FormiJation  of  BEM  equations  requires  numerous  Green's  function 
evaluations,  and  thus  direct  numerical  quadrature  is  simply  too 
computationally  expensive. 

Previously,  the  best  available  algorithm  has  been  the  method  proposed  by 
Wilson  and  Cruse,  [1978].  In  their  approach,  reasonable  computation  time 
was  obtained  by  employing  cubic  interpolation  with  tabulated 
precalculated  values.  This  method  has  three  potential  shortcomings. 

First,  storage  requirements  can  become  significant,  especially  if  first  and 
second  derivatives  of  the  Green's  function  are  required  (such  as  for 
evaluation  of  hypersingular  equations  used  for  fracture  mechanics). 
Second,  as  noted  in  Scholar  [1994],  accuracy  may  be  a  concern,  as  the 
variation  in  the  function  values  for  a  highly  anisotropic  material  is  likely 
to  be  significantly  larger  than  for  a  mildly  anisotropic  solid,  and  thus  the 
number  of  precalulated  points  needed  to  maintain  accuracy  is  a  function  of 
the  anisotropy  of  the  material.  Third,  the  precalculated  points  are 
material  specific.  The  precalculation  process  must  be  performed,  and  the 
resulting  tables  must  be  maintained,  for  each  material  of  interest. 
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For  this  project  we  decided  to  implement  the  residue  calculus  technique 
for  evaluating  the  function  described  by  Sales  and  Gray  [1996].  This 
approach  gives  us  the  obvious  benefit  of  the  direct  participation  of  one  of 
the  authors,  access  to  the  software  developed  for  the  paper,  and  capability 
to  develop  a  subroutine  that  would  work  for  all  orthotropic  materi^s  with 
no  precedculations  or  storage  of  tables  required.  The  residue  approach  has 
also  been  exploited  by  Wang  [1995];  however,  this  paper  is  primarily 
theoretical,  and  does  not  focus  on  computations. 

Following  Sales  and  Gray,  the  contour  integral  above  is  written  as 

=  (3-7) 


where  ^  is  written  parametrically  as 


■|,W 

sin(6)cos(O  +  cos(0)cos(v/^)sm(O 

«<)= 

= 

-cos(6)cos(0  +  sin($)cos(y)sin(t) 

A(‘). 

-sin(y/)sin(t) 

The  integration  can  be  performed  using  residue  calculus  by  transforming 
the  integrand  into  a  rational  function  and  expanding  the  range  of  the 
integral  to  include  all  real  numbers.  As  shown  by  Dederichs  and  LeibfHed 
[1969],  this  can  be  done  with  the  substitution 

Z  =  tan(0.  (3-9) 

With  this  substitution,  ^  can  be  written  as 

^  =  cos(05(Z),  (3-10) 


where 


sin(0)  +  Z  cos(6)  cos(v'^) 


5(Z)  = 


-  cos(0)  +  Z  sin(9)  cosiy) 


-Zsin(vA) 


(3-11) 


Since  the  elements  of  ^  are  linear  functions  of  Z,  and  A  is  a  quadratic 
function  of  the  elements  of  K  are  second  order  polynomials  in  Z.  Thus, 
by  Cramer's  rule,  the  elements  of  K'^  are  rational  functions  of  Z, 


cos\t)  Q(Z)  ’ 


(3-12) 


where  P  and  Q  are  polynomials  of  Z, 
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=  (3-13) 

/=0 

Q(Z)  =  jlb,Z‘.  (3-14) 

1=0 

With  the  observation  that  l/cos^(^)  =  dZjdt,  the  Green's  function  can  be 
written  in  the  form 


r~  P  AZ) 

=  l\  ^dZ 

(3-15) 

J-  Q(Z) 

=  4my  Residue ,  — - 

S  G(Z) 

(3-16) 

=  44W-A.) 

s  e(z) 

(3-17) 

ti  Q(K) 

(3-18) 

where  Xj,  Xj,  and  X3  represent  the  three  roots  of  Q(Z)  in  the  upper  half  of 
the  complex  plane. 

To  evaluate  the  Green's  function  numerically  using  this  technique  is  a 
two  step  process.  First,  the  three  roots,  X„,  are  determined  using  Newton's 
method.  Second,  Eq.  3-18  is  evaluated. 

The  coefScients  of  P  and  Q  and  6,  in  3-13  and  3-14)  are  functions  of  the 
components  of  K'S  and  thus  depend  on  the  material  constants.  For 
orthotropic  materials  there  are  nine  independent  material  constants 
(Cii,Ci2, ...  Cgg).  For  the  current  project,  the  Maple  V  symbolic  algebra 
program  was  used  to  construct  the  K  matrix  symbolically  in  terms  of  the 
nine  material  constants.  The  matrix  was  then  symbolically  inverted,  and 
expressions  for  the  coefficients  of  P  and  Q  were  developed.  While  these 
expressions  are  rather  tedious,  the  Maple  program  is  capable  of 
automatically  generating  FORTRAN  code  for  their  evaluation. 

A  similar  technique  can  be  used  to  determine  first  and  second  derivatives 
of  the  Green's  function.  This  process  is  detailed  in  the  paper  by  Sales  and 
Gray,  and  will  not  be  repeated  here. 


3.2.2  Testing 

We  have  performed  a  series  of  tests  to  assess  the  correctness  and  accuracy 
of  the  new  software.  For  example,  consider  the  evaluation  of  the 
normalized  Green's  function  along  vectors  lying  in  the  a  plane  passing 
through  the  y  axis  and  making  a  45°  angle  with  the  x  and  z  axis,  as  shown 
in  Figure  3-3. 
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We  have  evaluated  the  Green's  function  using  both  our  software,  and  a 
direct  numerical  integration  using  Simpson's  rule  (for  60  segments).  A 
plot  of  the  relative  percentage  error  for  two  of  the  components  is  shown  in 
Figure  3-4. 


Figure  3-3.  Path  for  the  evaluation  of  Green's  function 

As  c£m  be  seen  by  this  graph,  the  accuracy  of  this  approach  is  very  good, 
with  maximum  errors  less  than  1/lOOth  of  a  percent. 


Figure  3-4.  Relative  percent  error  for  two  components  of 
the  Green's  function 
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3^  Development  of  a  3D  Symmetric  Galerkin  BEM  Code 

A  substantial  portion  of  the  resources  for  this  project  were  spent 
developing  a  new,  3D,  symmetric-Galerkin  code  and  interfacing  it  with  our 
FRANC3D  modeling/meshing/postprocessing  system.  This  work  is 
forward-looking  towards  Phase  II,  where  it  is  anticipated  that  the 
symmetric-Galerkin  formulation  will  provide  an  efficient  boundary 
element  technique  that  can  accurately  determine  stresses  and 
displacements  near  cracks  and  other  boundaries. 


3.3.1  Code  Organization  and  Development 

There  are  three  major  tasks  that  a  BEM  program  is  required  to  perform: 
read  input  data  and  precompute  heavily  reused  data,  perform  all  element 
integrations,  and  solve  the  resulting  system  of  simultaneous  equations. 
The  software  in  our  new  program  is  organized  along  these  lines. 


Read  Input  and  Precompute  Element  Data 

The  subroutines  that  read  the  input  file  and  precompute  element 
information  are  written  in  the  C  programming  language.  We  feel  that  C 
language's  capacity  for  handling  dynamic  memory  allocation  and  for 
manipulating  character  based  information  makes  it  an  ideal  choice  for 
these  routines.  In  particular,  because  of  the  dynamic  memory 
management  the  language  provides,  there  is  no  need  to  "preset"  limits  on 
the  niunber  of  nodes,  elements,  and  materials  that  the  program  can 
handle.  Such  limits  will  be  determined  at  runtime  based  on  the  available 
resources. 

With  a  view  towards  flexibility  and  computational  efficiency,  there  are  two 
concepts  that  we  exploit  quite  heavily  in  the  code:  generalized  element 
data,  and  precomputation  of  this  data.  By  generalized  element  data  we 
mean  that  much  of  the  information  that  one  needs  on  an  element  basis  has 
been  formulated  in  a  manner  that  will  work  for  all  linear  and  quadratic 
order  elements.  For  example,  one  often  needs  to  determine  the  Cartesian 
coordinates  for  a  point  in  an  element  given  in  parametric  coordinates. 

This  is  easily  and  conventionally  done  by  using  the  element  shape 
functions  to  interpolate  the  element's  nodal  coordinates.  That  is, 
symbolically. 


nnodes 

Xa  =  a  =  1..3.  (3-19) 

1=1 

Computationally,  this  procedure  becomes  more  involved.  One  must  first 
compute  the  values  of  the  shape  functions  at  the  specified  parametric 
coordinates.  These  will  be  different  and  of  different  number  for  each  type 
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of  element.  This  leads  to  case  analyses  illustrated  in  the  following 
FORTRAN  code; 

if  (nnodes  .eq.  3)  then 

call  T3_shape (u,v, shape) 
else  if  (nnodes  .eq.  4)  then 
call  Q4_shape (u,v, shape) 
else  if 

end  if 


do  j  =  1,3 

coord(j)  =  0.0 
do  i  =  1, nnodes 
this_node  = 

1  node_coords ( j , connectivity (i, this_elem) ) 

coord(j)  =  coord(j)  +  shape (i) *this_node 
end  do 
end  do 

Such  case  analyses  are  undesirable  for  two  reasons.  First,  if  a  new 
element  is  added  to  the  program  (e.g.,  a  new  crack  tip  element),  all  such 
case  analyses  must  be  found  and  modified.  Second,  such  case  analyses  can 
seriously  hinder  or  defeat  a  compiler’s  effort  to  optimize,  vectorize,  or 
parallelize  the  code. 

In  our  code,  this  evaluation  is  generalized  for  all  element  types  using  the 
expression, 

x„  =  +  a„2M  +  a„3V  +  a„y  +  a^^uv  +  +  a^^uv^,  a  =  1..3.  (3-20) 

The  corresponding  FORTRAN  code  is  thus, 

u2  =  u  *  u 
v2  =  V  *  V 
uv  =  u  *  V 
u2v  =  u2  *  V 
uv2  =  u  *  v2 


do  j  =  1,3 

coord(j)  =  a(j,l)  +  a(j,2)*u  +  a(j,3)*v  + 

1  a(j,4)*u2  +  a(j,5)*uv  +  a(j,6)*v2  + 

1  a(j,7)*u2v  +  a(j,8)*uv2 

end  do 

It  may  appear  that  this  computation  would  be  expensive  for  the  simpler 
elements,  where  many  of  the  a  coefficients  would  be  zero.  Actually,  if  one 
accounts  for  the  computation  of  the  shape  function,  and  the  loop  and 
conditional  overhead,  the  latter  is  only  slightly  more  expensive  in  terms  of 
operations  for  these  elements,  and  significantly  cheaper  for  the  higher 
order  elements. 
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This  example  is  a  simplification  of  the  actual  situation  in  the  code.  There 
we  need  to  be  able  to  evaluate  the  coordinates  for  point  in  terms  of  a  polar 
expansion  (p,0)  about  any  location  in  the  element.  The  resulting 
expression  is, 


P  + 


=  («al  +  dalU  +  ^<*3^  +  +  a„,uv^)  + 

K9  +  +  fl„,4V^)cOS(0)  + 

_(a„i5  +  +  fl„,7V  +  )sin(e)  _ 

(a„2,  +  a„^u  +  a„23v)cos^(0)  + 

(a„24  +  «a25“  +  a<,26v)cos(0)sin(0)  + 

L(^a27  +  «a28«  +  fla29v)sin^(0) 

[(a„3o)cos^(0)sin(0)  +  (a„3,)cos(P)sin^(0)]p^ 


(3-21) 


|P'  + 


There  are  two  classes  of  precomputed  data,  those  which  are  done  for  each 

individual  element,  and  those  which  are  done  for  each  element  type. 

The  following  data  is  precomputed  in  the  program: 

for  each  element 

1.  The  Cartesian  coordinates  of  all  Gauss  integration  points, 

2.  The  integration  weights  for  all  Gauss  points, 

3.  The  element  normal  for  all  Gauss  points, 

4.  The  coefficients  needed  to  find  the  Cartesian  coordinate  for  any  location 
in  an  element  given  the  parametric  coordinates  of  an  origin  point  and 
polar  offsets  from  that  origin  (Eq.  3-21),  and 

5.  The  coefficients  needed  to  evaluate  the  Jacobian  at  any  location  in  an 
element  given  the  parametric  coordinates  of  the  point. 


for  each  type  of  element 

1.  The  parametric  coordinates  of  the  Gauss  integration  points, 

2.  The  integration  weights  for  all  Gauss  points, 

3.  The  shape  function  values  at  all  Gauss  points, 

4.  The  beginning  £md  ending  angle  and  the  coefficients  needed  to  find  the 
distance  to  the  element  edge,  for  a  polar  integration  over  the  element 
area,  centered  at  each  comer  node, 

5.  The  coefficients  needed  to  find  the  values  of  the  shape  functions  at  any 
location  in  the  element  given  the  parametric  coordinates  of  an  origin 
point  and  polar  offsets  from  that  origin,  and 

6.  The  parametric  coordinates  of  the  comer  nodes  of  the  element. 
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By  precomputing  all  this  information,  we  have  completely  avoided  any 
case  analysis  inside  any  of  the  element  integration  routines. 

Furthermore,  since  much  of  this  information  is  heavily  reused  during  the 
integration  process,  we  save  the  cost  of  recomputing  the  information  every 
time  it  is  needed. 


Element  Integrations 

Much  of  the  program  is  devoted  to  computing  element  integrations.  This 
portion  of  the  program  is  written  in  FORTRAN.  FORTRAN  is  the 
conventional  programming  language  for  such  computationally  intense 
algorithms,  and,  in  general,  the  best  optimizing  and  parallelizing 
compilers  available  are  for  FORTRAN  code. 

There  are  four  different  cases  of  element  integrations  depending  on  the 
physical  locations  of  the  source  and  destination  elements.  These  are 
illustrated  in  Figure  3-5.  In  all  cases  the  integrands  involve  1/r  or  l/r® 
terms,  leading  to  singular  integrals  in  three  of  the  four  cases. 

Most  integrations  are  of  the  Disjoint  type  (Figure  3-5a).  This  is  a 
nonsingular  integration.  Coincident  integrations  (Figure  3-5b)  occur 
when  the  source  and  destination  element  are  the  same.  This  integral 
becomes  singular  as  the  Q  point  approaches  the  P  point.  There  are  two 
possible  types  of  adjacent  integrations.  When  the  two  elements  share  a 
common  edge  (Figure  3-5c)  the  integrand  is  singular  along  the  entire 
common  edge.  W^en  the  two  elements  share  a  common  vertex  (Figure  3- 
5d)  the  integrand  is  singular  only  at  this  vertex. 

Within  the  code,  different  computational  strategies  are  employed  for  each 
of  the  four  cases.  The  approaches  are  summarized  here  and  described  in 
detail  in  Appendix  A. 


Disjoint 

This  case  is  relatively  straightforward  and  is  handled  with  conventional 
Gauss  integration.  The  integrand  contains  products  of  the  element  shape 
function  values,  Jacobean  determinates.  Gauss  integration  weights,  and 
appropriate  kernel  functions. 


Coincident 

The  integrand  in  this  case  becomes  singular  as  the  source  point 
approaches  the  destination  point.  To  handle  this  situation,  the  integrand 
is  decomposed  into  singular  and  nonsingular  parts.  If  for  each  gauss  point 
a  transformation  is  made  to  a  polar  coordinate  system  centered  at  the 
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Gauss  point,  an  analyticeJ  expression  for  the  singular  part  of  the 
integrand  can  be  found. 


Figure  3-5  The  four  possible  configurations  for  element  integrations. 


The  first  step  is  to  break  the  element  into  a  number  of  triangles,  4  as 

shown  in  Figure  3-6  (the  description  given  here  is  for  quadrilateral 
elements,  but  works  similarly  for  triangular  elements).  Within  each 
triangle,  the  integration  over  6  is  performed  numerically  (that  is,  discrete 
angles  and  weights  are  chosen  based  on  standard  quadrature  rules).  The 
integration  over  p,  however  is  performed  analytically.  This  process  is 
repeated  for  each  of  the  Gauss  integration  points  in  the  element. 


Vertex  Adjacent 

The  integrand  in  this  case  becomes  singular  at  the  common  vertex.  To 
handle  this  situation,  the  integrand  is  again  decomposed  into  singuleir  and 
nonsingular  parts.  Again,  an  analytical  expression  for  the  singular 
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integral  can  be  found  for  a  polar  coordinate  system,  in  this  case,  centered 
at  the  common  vertex. 


Figure  3-6  Coincident  polar  element  integration. 


This  situation  is  illustrated  in  Figure  3-7.  The  Q  element  is  decomposed 
into  two  triangles,  and  Tj  (only  one  triangle  is  necessary  for  triangular 
elements).  Within  each  triangle,  the  integration  over  0^  is  performed 
numerically.  The  integration  over  p  ,  is  partly  analytic  and  partly 
numerical. 


Figure  3-7  Vertex  adjacent  element  integration. 

A  polar  coordinate  system  centered  at  the  common  vertex  is  used  in  the  P 
element  also.  However,  in  this  case,  both  the  and  integrations  are 
done  numerically  by  using  the  standard  Gauss  integration  points  emd 


31 


F33615-96-C-5070  SBIR  Topic  AF96-150:  Phase  I 


Final  RepoH 


weights  (and  finding  the  appropriate  0^  and  that  correspond  to  these 
points). 


Edge  Adjacent 

The  edge  adjacent  case  is  similar  to  vertex  adjacent,  but  now  the 
singularity  is  along  the  entire  common  edge.  The  integrand  is  again 
decomposed  into  singular  and  nonsingular  parts.  In  this  case,  however, 
we  use  mixed  coordinate  systems. 

As  illustrated  in  Figure  3-8,  a  polar  coordinate  system  is  used  for  the  Q 
element.  The  origin  of  the  coordinate  system  is  a  point  on  the  common 
edge.  This  point  is  the  projection  of  the  P  integration  point  onto  this  edge, 
and  thus  moves  for  the  different  Gauss  points  in  the  P  integration. 


Figxire  3-8.  Edge  adjacent  element  integration. 

The  Q  element  is  decomposed  into  triangles,  two  for  triangular  element 
and  three  for  quadrilaterals.  Within  each  triangle,  the  integration  over  tf 
is  performed  numerically.  The  integration  over  p^,  is  performed 
analytically. 

The  P  element  integration  is  performed  in  a  rectilinear  coordinate  system. 
The  integration  for  the  u  coordinate  is  performed  numerically;  however,  an 
anal3rtical  expression  has  been  derived  for  the  integration  along  the  v 
coordinate. 
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Our  Eilgorithm  for  treating  the  edge  adjacent  case  is  a  new  development. 
The  resxilting  analytical  expressions  are  complicated,  but  £u*e  easily 
derived  using  S3mibolic  computations. 


Solving  Systems  of  Equations 

Once  element  integrations  have  been  performed,  the  resulting  system  of 
simultaneous  linear  equations  must  be  solved.  While  this  procedure  is 
computationally  expensive,  it  is  also  well  imderstood,  and  many  linear 
equation  solves  are  available.  Our  intention  is  to  not  spend  any  more 
effort  on  this  part  of  the  program  than  necessary. 

For  testing,  where  the  coefficient  matrix  is  small  enough  to  be  contained  in 
virtual  memory  without  excessive  thrashing,  we  are  using  the  widely 
available  LINPAC  solver. 


When  we  start  working  with  realistic  problems,  with  correspondingly 
larger  coefficient  matrices,  we  anticipate  using  the  out-of-core  BEM 
optimized  solver  that  we  developed  for  our  existing  BEM  program. 


3.3.2  Integration  with  FRANC3D 

We  have  modified  the  FRANC3D  program  to  generate  input  files  for  the 
new  BEM  code.  This  effort  was  more  involved  than  originally  anticipated 
because  the  new  code  treats  "corner"  boimdary  conditions  differently  than 
our  existing  BEM  code. 

In  the  new  program,  "ambiguous"  comer  boimdary  conditions  are  treated 
by  duplicating  the  appropriate  nodes.  For  each  degree-of-freedom  at  each 
node,  either  a  traction  or  displacement  can  be  specified,  with  the  other 
being  an  unknown.  At  comers,  however,  traction  values  are  not 
necessarily  unique,  leading  to  a  possibility  of  more  than  one  unknown  at 
the  comers.  If  this  is  the  case,  we  need  to  "duplicate"  the  node  so  that  the 
number  of  knowns  matches  the  number  of  unknowns.  Duplicating  the 
node  affects  the  node  list  and  the  nodal  connectivity  table,  thus 
complicating  input  file  generation. 

This  situation  can  be  illustrated  in  2D.  Figure  3-9  shows  (a)  a  comer  and 
(b)  a  non-comer  configuration.  First  assume  that  the  displacement  at  the 
node  is  known.  In  case  (b),  where  ni  equals  /i2,  ti  equals  t2,  we  have  one 
unknown  (traction)  for  our  known  (displacement).  In  case  (a),  however, 
the  normals  and  corresponding  tractions  will  be  different.  Therefore,  we 
have  two  unknown  (tractions)  for  one  known  (displacement).  In  this  case 
we  must  duplicate  the  node  so  that  both  tractions  appear  as  unknowns  in 
the  system. 
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Figure  3-9.  2D  comer  and  non-corner  boundary  condition  configurations. 

If,  on  the  other  hand,  tractions  are  known  at  the  comer,  there  are  again 
two  possible  cases.  If  the  tractions  are  equal  then  we  have  one  known  and 
one  unknown.  If,  however,  the  tractions  are  not  equal,  there  are  two 
knowns  for  one  unknown,  and  we  must  duplicate  the  node. 

The  3D  case  is  a  generalization  of  this  2D  analysis,  but  is  significantly 
complicated  by  the  fact  that  there  may  be  any  number  of  elements  incident 
at  a  node,  each  with  their  own  known  and  unknown  boundary  condition 
type  and  value.  In  the  3D  case,  the  actual  number  of  nodes  created  will  be 
between  one  and  the  niimber  of  incident  elements. 

To  determine  the  number  of  actual  nodes  created  in  the  input  file  for  any 
node  location,  the  following  general  algorithm  was  developed: 

1.  Assume  initially  the  worst  case  that  one  node  will  be  generated  for  each 
incident  element.  We  generate  a  data  stmcture  for  each  such  "use"  of  a 
node,  which  contains  the  types  of  known  and  unknown  boundary 
conditions,  the  known  boundary  condition  values,  and  the  element's 
normal. 

2.  Compare  all  possible  combinations  of  node  uses  and  determine  which 
node  uses  have  equivalent  knowns  and  unknowns,  and  can  be  grouped 
together. 

3.  Two  node  uses  are  determined  to  be  equivalent  if  they  meet  the 
following  criteria: 

a.  They  both  have  the  same  types  of  known  boundary  conditions  for 
all  three  coordinate  degrees  of  freedom. 

b.  If  the  tractions  are  known,  then  all  three  traction  values  must  be 
equal. 

c.  If  the  displacement  is  known,  then  the  element  normals  must  be 
equal. 

4.  One  node  is  generated  in  the  input  file  for  each  group  of  element  uses, 
emd  the  appropriate  element  nodal  connectivity  is  generated. 
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The  above  algorithm  has  been  implemented  and  tested  within  FRANC3D, 
and  appears  to  be  working  well. 

The  remainder  of  the  BEM  input  file  creation  is  a  relatively  streught 
forward  specification  of  nodal  coordinates,  element  connectivity,  material 
properties,  and  boimdary  conditions.  The  only  significant  additional 
consideration  is  the  case  of  bi-material  interfaces.  For  the  elements  on  the 
interface,  the  BEM  code  expects  to  see  two  elements  in  the  input  file,  with 
reversed  nodal  ordering  (each  ordering  corresponds  to  an  outward 
pointing  normal  for  the  two  materials)  and  differing  material  numbers. 
This  is  easily  handled  by  way  of  the  FRANC3D  internal  data  structures. 


4.  Summary  and  Conclusions 

We  have  successfully  completed  a  six  month  Phase  I  project  to  show  the 
feasibility  of  developing  a  3D  boundary  element  analysis  system  for 
composite  joints  with  discrete  damage.  Our  overall  goal  for  Phase  I  was  to 
demonstrate  expertise  in  BEM  methods  appropriate  to  the  detailed  stress 
analysis  of  arbitrary  bolted  or  bonded  composite  joints  containing  cracks 
and  contact  surfaces. 

We  have  completed  all  three  tasks  in  the  Phase  I  workstatement: 

•  We  augmented  the  Galerkin-based,  2D,  orthotropic  BEM  program 
developed  for  the  original  program  solicitation  to  include  hypersingular 
BEM  traction  equations  so  that  a  S3nmmetric  coefficient  matrix  could  be 
formxilated.  We  believe  that  this  orthotropic  symmetric-Galerkin  code 
is  imique, 

•  Using  the  code  developed  in  Task  1,  we  solved  an  example  laminate 
problem  [Pagano  1978],  to  show  that  the  symmetric-Galerkin  (SG) 
approach  can  yield  accurate  interlaminar  stress  predictions  near  free 
surfaces. 

•  We  demonstrated  a  capability  to  perform  3D  contact  analysis  using  our 
existing  3D  BEM  program.  The  approach  is  a  penalty  function 
technique  derived  from  a  current  capabilities  for  fluid  driven  fracture, 
and  can  also  be  used  for  the  proposed  Phase  II  3D  simulation  system. 

In  addition  to  these  tasks,  considerable  resources  were  spent  on  three 
additional  tasks  that  further  our  goal  of  demonstrating  expertise,  and  are 
preparatory  for  a  potential  Phase  II  effort. 

•  We  have  added  fracture  analysis  capabilities  to  the  2D  SG  code  of  Task 
1.  We  believe  that  this  is  the  first  SG  fracture  for  elasto-statics,  and 
certainly  the  first  for  orthotropic  materials. 
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•  We  have  developed  software  to  evaluate  the  3D  anisotropic  Green's 
function.  This  is  a  generalization  of  the  work  of  Sales  £ind  Gray  [1996]. 
In  this  approach  residue  calculus  is  used  to  transform  the  Green's 
function  to  a  ratio  of  pol3momials.  A  symbolic  algebra  program  is  used 
to  generate  polynomial  coefficients.  Newton's  method  is  used  to  find  the 
roots  of  the  polynomial.  Our  implementation  can  handle  multiple 
orthotropic  materials;  no  precomputations  are  required. 

•  We  have  developed  a  fully  3D  Galerkin  BEM  program.  Coding  for  this 
program  has  been  completed,  we  are  currently  in  the  testing  and 
debugging  phase.  This  code  includes  linear  and  quadratic  order 
elements,  and  has  been  interfaced  with  our  FRANC3D  program  that  we 
use  to  generate  models,  meshes,  and  boimdary  conditions. 

In  summary,  we  have  successfiilly  completed  the  work  tasks  specified  in 
our  proposal.  In  addition,  we  have  gone  substantially  beyond  the  original 
Phase  I  objectives  towards  creating  a  capability  for  3D  boundary  element 
analysis  of  composite  joints  with  discrete  damage.  We  believe  that  we 
have  demonstrated  that  the  sjnnmetric-Galerkin  BEM  formulation  is  well 
suited  for  this  task,  and  that  Fracture  Analysis  Consultants  has  the 
expertise  to  pursue  this  objective. 
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Appendix  A  -  Details  of  the  3D  SG-BEM  Formulation 

Symbols: 

Uj  -  coefficients  of  linear  term  of  Taylor  expansion  of  distance 
between  p  and  q,  (j  =  1..3) 
a  -  RMS  norm  of  the  a  vector 
C  -  material  constitutive  matrix  (6  x  6) 

-  matrix  of  derivatives  of  the  traction  kernel  (kk  =  1..3,  jj  =  1..3, 
ll  =  1..3) 

-  matrix  of  derivatives  of  the  traction  kernel  for  a  truncated  Taylor 
expanion  in  polar  coordinates  (kk  =  1..3,  jj  =  1..3,  ll  =  1..3) 

Ep,Eq  -  elements  for  integration 

-  normalized  Green’s  function  (kk  =  =  1..3) 

Gu^jj  -  truncated  Taylor  expansion  of 

gp  -  Gauss  point  in  P  element 

gq  -  Gauss  point  in  Q  element 

|/|  -  Jacobian  determinant 

PI  -  truncated  Taylor  expansion  of  |/|  in  polar  coordinates 
L  -  distance  from  gauss  point  to  element  edge  for  element 
singxilar  polar  integrations 
P  -  integration  source  point 

p^  -  integration  source  point  Cartesian  coordinates  (jj  =  1..3) 

Q  -  integration  field  point 

Qjj  -  integration  field  point  Cartesian  coordinates  (jj  =  1..3) 

Qjj  -  field  point  coordintates  for  a  tnmcated  Taylor  expansion  of  r  in 
polar  coordinates  (jj  =  1..3) 

r  -  distance  between  source  and  destination  point 

f  -  truncated  Taylor  expansion  of  r  in  polar  coordinates 

'^kkjj  ■  displacement  kernel  (kk  -  \..3,jj  =  1..3) 
fu^jj  -  truncated  Taylor  expansion  of  Ttkjj  in  polar  coordinates 
Uu,jj  -  traction  kernel  (kk  =  \..3,jj  =  1..3) 

-  displacement  vector  (kk  =  1..3) 

w  -  Gauss  integration  weight 


41 


F3.1615-96-C-5070  SBTR  Topic  AF96-150:  Phase  I 


Final  HepoH 


p-  -  polynomial  coefficients  for  a  polar  expansion  of  the  shape 
functions  (i  =  0..3) 

-  vector  of  strain  vectors  Qj  =  1..3) 

ijj  -  vector  of  strain  vectors  for  truncated  Taylor  expansion  in  polar 
coordinates  (jj  =  1..3) 

Yi  -  polynomial  coefficients  for  a  polar  expansion  of  r  (i  =  0..4) 

Gjj  -  vector  of  stress  vectors  (jj  =  1..3) 

Ujj  -  vector  of  stress  vectors  for  truncated  Taylor  expansion  in  polar 

coordinates  (JJ  =  1..3) 

-  traction  vector  (kk  =  1..3) 

-  shape  function  (k  =  1.. number  of  element  nodes) 

9  -  angle  for  polar  and  sherical  coordinates 

^  -  angle  for  spherical  coordinates 

p  -  distance  for  polar  coordinates 

77,4  -  p  element  parametric  coordinates 

-  q  element  parametric  coordinates 


Governing  Equations: 

The  governing  boimdary  element  equation  is 

+  J u,(Q)  dQ  =  J U^,t„(Q)dQ  . 

In  terms  of  a  Galerkin  formulation,  this  can  be  rewritten  as 

I V'. (/>)«„(/>)  +  J V,(P)\ T^ji(F.Q) dQ  =  j r,(P)f  U^j,Tj,{Q)dQ  , 

with  the  traction  kernel  function  defined,  in  terms  of  spherical 
coordinates,  as 


(1) 

(2) 

(3) 


The  displacement  kernel  fimction,  T,  is  derived  from  the  derivative  of  U 
with  respect  to  the  three  coordinate  directions. 


dU. 


kkjj.ii 


%K 


Qii-Pu 

-3 


P,  ,  1  ^kkjj  d6  idGi 
r  dd  dq„ 


_ kkjj  dip 

r  d<f)  dq„ 


(4) 
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from  which  can  be  formulated  a  strain  vector 


£y2  1  (^^2i;,3  ) 

Using  the  constitutive  relationship 


dU^JJ,^ 

dU.jj,, 


^jj  -  Qx6  ^jj 


the  displacement  kernel  can  be  defined, 


- 

Tkkjj  =  (^xy 

Spherical  coordinates  are  usually  defined  as  follows: 


-Pi  (I2-P2  ^3-P3  1- 


r<?i-A 

I  r  ’ 


=  (cos  d  sin  <t),  sin  d  sin  0,  cos  0) 


V  r  r  rj 

with  O<0<;r,  O<0<27r.  In  this  case, 


de  _  q2-P2  1 


de  _Qi-Pi  1  i^_0  and 


r  sirr  (ft  dq^  r  sin  <j>  dq^ 


^  i^. =k:;£j)f2!i,i^.=-f5L^coseH-^:^sine\  ao) 

dq^  r  sin0  dq^  r  sin^  dq.^  \  r  r  ) 

However,  this  gives  rise  to  a  singular  point  along  the  z  axis.  Therefore,  in 
cases  where  {q^  -  p^lr  is  >  0.7071,  the  following  alternative  spherical 
coordinate  system  defintion  is  used, 


— h.  Bl — El.  Si — =  (cos0,cos0sin^,sin0sin 0), 
for  which. 
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_Qi-P2  1 

2  •  1  2  •2j>  cUlU 

dq^  r  sin  ^  dq^  r  sin  0 


(12) 


^  =  -r«^acose  +  ^sin9\-^  =  ^^5!£,  i^ij^cose 
V  ^  r  J  dq2  r  sin0  ^3  r  sin0 


(13) 


Single  integral  term 

The  single  integral  term  can  be  evaluated  numerically  as  follows, 
j  Wk(P)^kk^Pj)dP  =  u^'^\if^{gp)\ifj{gp)\J{gp\w{gp).  (14) 

gp 

Nonsingular  integrals 

Nonsingular  integrals  arise  when  elements  P  and  Q  have  no  nodes  in 
common.  In  this  case  both  integrals  can  be  evaluated  numerically;  the 
traction  integral  being: 

l,Wk(P)lu^jj(P,Q)dP 

=  vMdnd^  Wj  Ui^Md^  dl 

=  Yy^kigPppigpiwigpy^  U^jjigP.gQ)  YjigQpoigQpigq) 

gp  8<l 

= ^'Ly^i^^sppp^spiwigpyZ  wM(ipQ^g<i)\w(<g(i) 

gp  gq  ^ 

and  the  displacement  integral  is 

=  ^^^Ly^i^^gPpp^SPiwpigP)  ^ftkjiigP^g(!)YjigQ)wigQ) 

gp  gq 

Note  that  in  this  case,  the  Q  element  Jacobian  determinant  is  included  in 
the  normal  vector  used  to  evaluate  the  displacement  kernel,  T  Eq.  (7). 


(15) 


(16) 


Coincident  Integrations 

Coincident  integrations  occur  when  P  and  Q  are  the  same  element.  In  this 
case  the  integrations  are  singular,  so  they  are  decomposed  into  a 
nonsingular  portion,  solved  numerically,  and  a  singular  portion,  solved 
analytically.  In  order  to  find  a  form  of  the  singular  integral  that  can  be 
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evaluated  anal3dically,  the  normal  area  integration  is  transformed  to  a 
polar  integration  centered  at  each  of  the  elements  gauss  points. 


To  do  this,  define  the  Q  gauss  points  as  polar  expemsions  about  the  P 
points,  which  are  located  at  the  element's  gauss  integration  points. 


jj  =  r}  +  pcosO,  and  ^=^  +  psind.  (17) 

One  can  then  define  the  vector  from  p  to  p  in  terms  of  p, 

^a-Pa=«aP  +  ^aP^+^^aP^  « =  U,3  (18) 

=  (of  +a\+  al)p^  +  2{a^by  +  +  a^b^)p^  + 

(2fl,c,  +  bi  +  la^c^  +  +  2a^c^  +  bl)p*  +  (19) 

2(feiC,  +  fejCj  +  b^c^)p^  +  icl  +cl  +  cl)p^ 

or 

=  /oP"  +  /iP'  +  Y2P'  +  YiP'  +  YaP'  •  (20) 

If  we  look  at  the  first  term  in  this  expansion,  o  =  +  d^-¥  a] ,  one  can  define 

^a=Pa+«aP>  (21) 

f  =  ap, 

(22) 

|/g  =|/g|^^o  “  k/’(^»^)|>  (23) 


Pu-Pii^ 
p3  ^kkjj 


A  A 

^  1  ^kkjj  do  ^  1  ^kkjj 
r  do  dq„  r  d<p  dq„ 


(24) 

(25) 
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a  =  Ce,  and 


(27) 


f  =  %h.  (28) 

Then,  the  traction  integral  becomes, 

=  r.k/.M'l  Vi  \  C„,|4|  dri  dl 

=I^Sv'*(«p)i'/p(sp)|'*’Csp)Jj^r;^6„|4|pdpde  (29) 

1  Pi  1  ^  I  |** 

=  ■^'L^kiSP)\JpigP)\^i8P)\^^ \¥j -Gktjj\jQ\ - ¥j jG^aYol  pdpdO 

The  first  part  is  well  behaved,  and  can  be  integrated  numerically, 

»yr  jp  „  Lr  r  j 

and  the  second  part  is  solved  with  a  hybrid  numeric/analytical  approach. 
Expanding  the  shape  functions  in  terms  of  polar  coordinates, 

=  )3o + A(0)P + P2ie)p^ + P3(e)p^  (31) 

then, 

-^^¥k(gP)\Jp(8Pi^(8P)  jldej^'^¥j^G^j\jQ\  pdp 
6K  gp  r 

=  ■^\h\L¥M\Jpi8P)H8P)  C^ddj^^'^Yj  pdp 

(32) 

=  -^\jQ\LyM\Jp(8P)\^(8P)  X 

r  a  L  2  3  4  J 

The  displacement  part  is: 
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(33) 


Y,(P)dPl^  Wj(Q)T^jj(PM)  dQ 

=  -^'^VMP)\Jpi8P)\^igP)\^^  V^7tti,ke|  PdpdB 

=  -^Ya^M\jpi8P)\^igP)\^  [v^y^«i/kG|  -  V^/tti/j4|]  P^P^^ 

oTC  gp  a 

+  -^^¥k(8P)\Jp(gP)HsP)j^^  wAampdpde 

Again,  the  first  part  is  well  behaved,  and  can  be  integrated  numerically, 

S  ■  Tu,jPq^.  (34) 

gP  8^ 

The  second  part  is  solved  with  a  hybrid  numeric/analytical  approach, 

-^J,¥M\jp(gP)H8P)  P^P 

gp 

=  ■^\jQ\L¥g(gP)\Jp(gP'i^(gP)  [dej^^'^fgg^Yjpdp  (35) 

gp 

=  -^m'L¥g(gP)\Jp(gPi^(gP)'Z^(^^)  C^^PiucjjYj  pdp 

gp  Afi 

rL(0)  ^ 

where  T^^jjYjpdp  is  evaluated  using  symbolic  computation. 

Vertex  Adjacent  Integrations 

When  the  two  elements  share  one  common  node,  the  integration  is  said  to 
be  "vertex  adjacent".  For  this  case,  the  traction  integral  is: 

^i¥,(P)dplYj(Q)l6^,(P,Q)<iQ 

=  -  Vi-Ai,]lQ  (36) 

The  first,  well  behaved,  part  is  integrated  numerically. 
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V't  ^gQ)~GitjjigP,gq)\j,  (^^)|  yv{gq) 

\  .  ”  1  -  .  (37) 

“  igP,gq)\j,(.gq)\^{gq) 

gp  gq  ^ 

The  second  part  is  evaluated  with  a  hybrid  numerical/analytical  approach, 


1  f®f2 


jj"  \l'^^^\j(<PpA)Ppdppddp\'^2 Vt(Pe.^e)^pG  ^Pg 


=  C  V;(P/.>^/>)P/>dp,.dV^p£'“^*“  V*(Pg.^o)^  ^Po 

=  -^^^^g^p'ppig^p'^<g^Q)\jQig^Q^^igPp)¥MPp^g^p)pP  X 

A 

|,'’“'"V.(Pe.9e)%-Pe‘^Pe 


with  the  analytical  expressions  for£^^^®Vt(Pe»^G)“^Po^PG 

determined  using  symbolic  computations,  the  resulting  expressions  are 
too  long  to  be  repeated  here. 

For  the  displacement  integral, 

the  first  part  is  integrated  numerically, 

i  X  v'*  (^p)  Vp  ^gpi  ^(gp)X  yj  (gp^gq)\Jq  (^7)|  w(gq) 

gp  gq  , 

1  ^  .  .  (' 

--^H¥i:igp)\jpigp)\^igp)'Z¥jigq)fujiigp^gq)\j,ig‘ii^^g^) 

8P  gq 

The  second  part  is: 
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V*(Pe’^e)  Pq  dpQ  dOQ 

^  'til  JeJ  V;(Pp.^/>)Pp^P/>^V^  pQ  ^PQ 


%Jt 


=  -^^w{g0p)\hi8^p'^^iS^Q)%^^Q^^^Pp')¥ MPp^g^p)Pp  X 

gOf  gOQ  gpr 


(41) 


gPr 
LgiOa) 


Jp  Wk  (Pq  ’  )  ^kkjj  Pq  dPa 


with  the  analytical  expressions  for  ¥k(PQ^^Q)'Pktjj Pq^Pq  being 
determined  using  symbolic  computations. 


Edge  Adjacent  Integrations 

The  traction  integral  for  this  case  is  the  same  as  Eq.  (xl).  Again  the  first 
part  is  evaluated  numericaly.  In  this  case,  however,  integration  over  the  u 
and  V  coordinates  in  the  P  element  are  done  independently. 


dQ 


=  ^  X  1'^/’ H  V'y  7  ke 


«« 


*V 


*? 


1  1  ^  I  ^  I 

“  E  ^(^“)E  ^gu,  gv)|7p  (gu,gv)\  wigv)^  Yj  (gq)  -  \jQ(gq)\  w(gq) 

07t  gu  gv  gq  ^ 


(42) 


The  analytical  part  in  this  case  is 


^  J  J  V^y(M,  v)dvdMj®'  l^^^Wk(P^^)jGkkji  pdpdd 

U  V 

=-^\Se\£''¥jMY,{p,e)jh^„pdpdvdddv 

T  '  .  . 

gB 

X  ^kkjj  ¥j  igu,  gv)  Yk  (P»  g&)^dp 

gv  ^ 


where  *V;(^w>^v)v/^*(p,g0)^dp  is  integrated  symbolically. 

Jo  ^  f* 


The  displacement  part  is  the  same  as  Eq.  (x3).  the  numerical  part  being. 
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=  -^'^w(gu)^V,,(gu,gv)\jp(gu,gv)\w(gv)Y,¥j(^q)T^j^\jQigq^wigq)  (44) 

OTT  gu  gv  gq 

--^'^w(gu)^Yk(g^rgv)\jp(^u,gv)\wigv)Y,Vj(gQ)faj^\jQigq^w(gq) 

OTC  gu  jv  gq 

The  analjiiical  part  is 

^  J  J  VjMdvduj'"  pdpdO 

uv  ' 

t  '  ,  . 

oTC  gu  gg 

X  fgg  jjy/j  (gu,  gv)  Yg(p,ge)  p  dp 

g'> 

with  j^^^^fujjYj(gu,gv)Yk(.p,gO)pdp  integrated  using  symbolic  computations. 
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